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The truncation of a pair potential at a distance r c is well-known to imply in general an impulsive 
correction to the pressure and other moments of the first derivatives of the potential. That depending 
on r c the truncation may also be of relevance to higher derivatives is shown theoretically for the Born 
contributions to the elastic moduli obtained using the stress-fluctuation formalism in d dimensions. 
Focusing on isotropic liquids for which the shear modulus G must vanish by construction, the 
predicted corrections are tested numerically for binary mixtures and polydisperse Lennard- Jones 
beads in, respectively, d — 3 and d — 2 dimensions. 
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I. INTRODUCTION 

Background. It is commonpractice in computational 
condensed matter physics [l|-[3| to truncate a pair inter- 
action potential U(r) at a conveniently chosen cutoff dis- 
tance r c with r being the distance between two particles 
i and j. This allows to reduce the number of interactions 
to be computed — energy or force calculations become 
thus (9(n)-processes with n being the particle number — 
but introduces some technical difficulties, e.g., instabili- 
ties in the numerical solution of differential equations as 
well-studied in the past especially for the molecular dy- 
namics (MD) method [j], |j] . Let us label the interaction 
between two beads i < j by an index I. For simplicity of 
the presentation and without restricting much in practice 
the generality of our results, it is assumed below that 

• the pair potential scales as U(r) = u(s) with the 
reduced dimensionless distance s — r/o~i where o~i 
characterizes the range of the interaction I and 

• the same reduced cutoff s c — r c /ai is set for all 
interactions I. 

For monodisperse beads with constant bead diameter a, 
as for the standard Lennard- Jones (LJ) potential [l[ 

u hJ (s) = 4e (Jg - 1) , (1) 

the scaling variable becomes simply s = rja and the 
reduced cutoff s c — r c /a. Even if the truncated potential 

ut(«) = u(s)H(s c - s) (2) 
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with H (s) being the Heaviside function [f| is taken by def- 
inition as the new system Hamiltonian, it is well known 
that impulsive corrections at the cutoff have to be taken 
into account in general for the pressure P and other mo- 
ments of the first derivatives of the potential [2J . These 
corrections can be avoided of course by considering a 
properly shifted potential Q 

u s (s) = {u{s)-u{s c ))H(s c -s) (3) 

as emphasized below (Sec. Ill A)) . 

Goal of presented work. In this report we wish to re- 
mind that the standard shifting of a truncated poten- 
tial is insufficient in general to avoid impulsive correc- 
tions for moments of second (and higher) derivatives of 
the potential. We demonstrate here that this is par- 
ticulary the case for the Born contribution C^f^ S (de- 
fined below) to the elastic moduli computed using the 
stress-fluctuation formalism described in great detail in 
the literature 0, 0413 • This should be of importance 
for the precise localization of the transition between dif- 
ferent thermodynamic phases using the elastic moduli, 
especially for liquid (G = 0) to solid (G > 0) transitions 
in network forming systems where the shear modulus G 
plays the role of an order parameter |l8j |. This is the 
case, e.g. for colloidal gels [19| , hyperbranched polymer 
chains with sticky end-groups [2(| or bridged networks 
of telechelic polymers in water-oil emulsions pll . [22| or 
living polymer- like micellar systems j23| . 

Outline. The paper is organized as follows: After re- 
minding first in Sec. Ill Al the known corrections for the 
pressure and similar first derivates of the potential, the 
impulsive correction for the general Born contribution 
C^ lS is stated in Sec. EH We describe then in Sec. EE] 
the corrections for the compression modulus K and the 
shear modulus G in isotropic systems. We comment on 
polydispersity effects and mixed potentials in Sec. Ill Dl 
Our results are rephrased in terms of the radial pair dis- 
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tribution function g(s) in Sec. Ill El which allows to pre- 
dict the asymptotic behavior for large s c . Section Hill 
gives some technical details on the two numerical model 
systems 0, HH| in d = 3 and d = 2 dimensions used to 
test our predictions in Sec. lIVI We consider in this paper 
the liquid high-temperature regime of both (in principle 
glass-forming) models where the shear modulus G must 
vanish [y, |26| , since this provides a simple reference for 
testing the predicted corrections. 

II. THEORETICAL PREDICTIONS 

A. Reminder 

Truncated potential. As usual for pairwise additive 
interactions the mean pressure P — Pd + P cx may be 
obtained as the sum of the ideal kinetic contribution 
Pd = k^Tp and the excess pressure contribution p5| 

with p — n/V being the number density, n the particle 
number, V the <i-dimensional volume and (. . .} indicating 
the usual thermal average over the configuration ensem- 
ble. (A prime denotes a derivative of a function with 
respect to its argument.) By taking the derivative of the 
truncated potential 

u[{s) = u'{s)H{s c - s) - u(s)5{s - s c ) (5) 

the excess pressure may be written as the sum P ex = 
P cx + AP 0X of an uncorrected (bare) contribution P cx 
and an impulsive correction AP CX . The latter correction 
is obtained numerically from Q 

AP CX = lim hi(s) with 

h i( s ) = ^Y^( s ' u ( s i) 6 ( s ' ~ s )) ( 6 ) 

being a weighted histogram. In practice, the proper limit 
s — > s~ may be replaced by setting s — s c for the his- 
togram. 

Shifted potential. The impulsive correction related to 
first derivatives of the truncated potential can be avoided 
by considering the shifted potential u s (s), Eq. ((3]), since 
u' s (s) = u'(s)H(s c — s). With this choice no impulsive 
correction arises either for similar observables such as, 
e.g., moments of the instantaneous excess pressure tensor 

pa/J = 1 V „» du s(si) 

= -^5><( Si )n>f. (7) 
l 

Here sf, ... stand for the spatial components of the re- 
duced distance between the particles, nf — sf /s;, . . . for 



the corresponding components of the normalized distance 
vector and Greek letters are used for the spatial coordi- 
nates a,/5,7, 6 = l,...,d. (We remind that the mean 
excess pressure P cx is the averaged trace over the instan- 
taneous excess pressure tensor, P cx =< Tr[P cx ' 9 ] > /d.) 
Specifically, if the potential is shifted, all impulsive cor- 
rections are avoided for the excess pressure fluctuations 

Cfis = _p v ((pgpg) - (p e f ) (j*f)) ( 8) 

(/3 = l/ksT being the inverse temperature) which give 
important contributions — especially for polymer-type 
liquids [IE U3] and amorphous solids [HI UH — to the 
elastic moduli computed using the stress-fluctuation for- 
malism 



B. Key point made 

Correction to the Born term. An even more impor- 
tant contribution to the elastic moduli (especially at high 
densities) is given by the Born term C^f lS already men- 
tioned in the Introduction [271 ]. Being a moment of the 
first and the second derivatives of the potential it is de- 
fined as @, m ee m 

CT & = ~ E ( (*?<(»0 - s t<M) nfnfnlnj) (9) 

using the notations given above. We remind that for 
solids with well-defined reference positions and displace- 
ment fields the Born contribution is known to describe 
the (free) energy change assuming an affine response to 
an imposed homogeneous strain [III 111,13 HE EH- As- 
suming now a truncated and shifted potential the impul- 
sive correction ACjf 7 * to C^ 5 = C^ S + ACf jS is 
simply obtained using 

<(s) = u"( S )H(s c - s) - u'(s)6(s - s c ) (10) 

which yields 

AC* P ~ fS = - lim hf l5 (s) with (11) 

s—>s c 

General impulsive correction. More generally, one 
might consider a property 

with f(s) being a specified function and (n) denoting 
the n-th derivative of the shifted potential u s (s). Let us 
further suppose that all potential derivatives up to the 
(n — 2)-th one do vanish at the cutoff s c . It thus follows 
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that A = A + AA takes an impulsive correction 
A A = — lim h n (s) with 

h n {s) = ^{/(oluMf,,)^,-,)) (12) 

being the relevant histogram. 

Generalized shifting. Obviously, the original poten- 
tial may not only be shifted by a constant u(s c ) but by 
a polynomial of s to make vanish the first and arbitrar- 
ily high derivatives of the potential at s — s c . In this 
way all impulsive corrections could be avoided in princi- 
ple. Since discontinuous forces at the cutoff may cause 
problems in MD simulations, a number of studies use for 
instance a "shifted-force potential" where a linear term is 
added to the potential [l|, |4| . The difference between the 
original potential and the generalized shifted potential re- 
moving the cutoff discontinuities means, of course, that 
the computed properties deviate to some extend from the 
original model. Only if the generalized shifting is weak, 
one may recover the correct thermodynamic properties 
using a first-order perturbation scheme pj. Since the 
(simply) shifted potential u s (s), Eq. (l3"l), is anyway the 
most common choice 0"S fl3 - fl7l . [24 . |25| , we restrict the 
presentation on this case and demonstrate how the im- 
pulsive correction associated to the non- vanishing u' s (s~) 
can readily be computed. 



C. Isotropic systems 

Lame coefficients. In order to show that these correc- 
tions may be of relevance we focus now on homogeneous 
and isotropic systems. We remind first that the two elas- 
tic Lame coefficients A and /i characterizing the elastic 
properties of such systems may be computed numerically 
using pj, 03 



prefactor stems from the assumed isotropy of the system 
and the mathematical formula 



A — Ap + Ab, 

M - -Pid = A*F + f-B 



(n>?y 



1 



d(d + 2) 



(1 + 26 a p) 



(15) 



(8 a f) being the Kronecker symbol pj) for the components 
of a unit vector in d dimensions pointing into arbitrary 
directions. Equation (|llj) implies then an impulsive cor- 
rection 

AAb = A/ie = — lim /i2(s) with 

h2{s) s d(dTW^^ u ' {si)s{si - s) ^ (16) 

Compression and shear modulus. Instead of using the 
Lame coefficients it is from the experimental point of 
view more natural to characterize isotropic bodies using 
the compression modulus K and the shear modulus G. 
The latter moduli may be expressed as 



K = (A + P) + -G, 

G = [1 - P = flB + - Pc: 



(17) 
(18) 



(13) 



We follow here the notation of Ref. [16] to emphasize the 
explicit pressure dependence which is often (incorrectly) 
omitted as clearly pointed out by Birch [29[ and Wallace 
[30l |. As one expects, kinetic elastic contributions terms 
do not enter explicitly for the shear modulus. Since only 
the Born contributions Ab = Mb cause a cutoff correction, 
this implies K = K + AK and G = G + AG with K and 
G being the uncorrected (bare) moduli and 

2 2 + d 

AK = AA B + 3 A MB - — — A^ B , (19) 
d d 

AG = A/i B (20) 

the impulsive corrections. We shall test these predictions 
numerically in Sec. IIVI 



where the only contribution due to the kinetic energy of 
the particles is contained by the ideal gas pressure Pid 
indicated for fi [28j. The first contributions indicated on 
the right hand-side of Eq. (TIB")) are the excess pressure 
fluctuation contributions Af and /if which may be ob- 
tained from the general C^ lS by setting, e.g., a — (3 — 1 
and 7 = S = 2 for Ap and a = 7 = 1 and (3 = 5 = 2 
for fip characterizing the shear stress fluctuations. The 
so-called "Born-Lame coefficients" flol 



A B = /x B = 



1 



d{d + 2)V 



£(( S ^"( Si )-W(sj))) (14) 



may be obtained from the general Born terms Cg' 37 ' 5 by 
setting, e.g., a = 7 = 1 and (3 = 6 = 2. The d-dependent 



D. Polydispersity and mixed potentials 

As stated in the Introduction we assume throughout 
this work the scaling U(r) = u(s) of the pair potential 
in terms of the reduced distance s = r/cri. This is done 
not only for dimensional reasons but, more importantly, 
to describe a broad range of model systems for mixtures 
and polydisperse systems where the interaction range 07 
may differ for each interaction I. Moreover, the type 
and/or the parameter set of the pair potential may vary 
for different interactions. For such mixed potentials u(s), 
Ut(s) and u s (s) and their derivatives take in principal an 
explicit index I, i.e. one should write u/(s), u t) ;(s), « s ,j(s) 
and so on. This is only not done here to keep a concise 
notation. For example one might wish to consider 
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• a generic polymer bead-spring model where some 
interactions I describe the bonded interactions be- 
tween monomers along the chain (which are nor- 
mally not truncated and need not to be corrected) 
and some the excluded volume interactions between 
the beads [17} . 

• the generalization of the monodisperse LJ poten- 
tial, Eq. (P, 

ui(s) = 4e; (s~ 12 - S - 6 ) with s = r/ai (21) 

where e; and o~i are fixed for each interaction I. In 
practise, each particle i may be characterized by 
an energy scale Ei and a "diameter" _Dj. The in- 
teraction parameters ei(Ei,Ej) and ai(Di,Dj) are 
then given in terms of specified functions of these 
properties. 

• the famous Kob- Andersen (KA) model for binary 
mixtures of beads of type A and B [24} , a partic- 
ular case of Eq. ([2~Tj) with fixed interaction ranges 
(TAAj 0bb and oab and energy parameters eaa, £bb 
and £ab characterizing, respectively, AA-, BB- and 
A£?-contacts. 

• a network forming emulsion of oil droplets in 
water bridged by telechelic polymers where the 
oil droplets are modeled as big LJ spheres, the 
telechelic polymers by a bead-spring model with a 
soluble "spacer" in the middle of the chain and in- 
soluble end-groups ( "stickers" ) strongly attracted 
by the oil droplets [2l|, |22| . Assuming sufficiently 
strong (in strength, number and life-time) sticker- 
oil interactions, such a system should behave as a 
soft solid with a finite shear modulus G (at least for 
a fixed finite sampling time) which may be probed, 
at least in principle, using Eq. (|18p. 



The impulsive corrections given in Eq. © for the pres- 
sure P, in Eq. (fTTj) for the general Born term C^f lS and 
in Eq. (p~6|) for the Born Lame coefficients Ab = Pb in 
isotropic systems stated all in terms of, respectively, the 
histograms hi(s), h^ 15 {s) and /12(a) remain indeed valid 
for such explicitly /-dependent potentials. From the nu- 
merical point of view, this is all what is needed and the di- 
rect computation of these histograms remains in all cases 
straightforward as illustrated in Sec. IIV AI 



E. Radial pair distribution function g(r) 

Notations. Especially for simple and c omp lex fluids 
and for all sorts of glass-forming systems 0, [26| it is com- 
mon practice to reexpress correlations and histograms in 
terms of the radial pair distribution function g{r) [J HJ. 
This is also of interest here since for large cutoff distances 
the pair distribution function must drop out, g(r c ) — » 1, 
allowing thus to predict the corrections in this limit. Let 
us remind first that, using the Gamma function T(x) |J, 



the (d — l)-dimensional surface of a d-sphere of radius r 
is given by 



Mr) 



r(d/2) 



for d = 2, 3, . . 



(22) 



and similarly for the (dimensionless) surface A(s) using 
the reduced distance s. 

Monodisperse interactions. For strictly monodisperse 
beads and similar interactions of constant interaction 
range a it is seen that Eq. ([6]) for the pressure correc- 
tion becomes 

AP CX = i ^p 2 a d A{s c )s c u(sc) x g(s c ) (23) 

where the factor 1/2 assures that every interaction is only 
counted once. Note that we have set g(s) = g{r/a) with- 
out introducing a new symbol. For the LJ potential, 
Eq. ((J), this leads to 



47T rf / 2 

T(d/2)d 



p 2 a d es d - 6 (l-s- 6 )xg(s c ). (24) 



Please note that the expression given in Ref. [2j is recov- 
ered by setting d — 3 and assuming g(s c ) ~ 1. Similarly, 
one obtains from Eq. (fT6|) the correction 



1 1 



2 d(d + 2) 



p 2 a d A(s c )siu'(s c ) x g(s c ) (25) 



for the Born-Lame coefficient we are mainly interested 
in. For a LJ potential this becomes 



Am 



24^ 2 



d(d + 2)T(d/2) 
where we have defined 



p 2 o- d ef hJ (s c ) x g(s c ) (26) 



/ LJ ( S )^(l-( So /s) 6 )/*' 



6-c 



(27) 



with so = 2 1 / 6 being the minimum of the potential. For 
sufficiently large cutoff distances where g(s c ) « 1 the 
correction thus decays as 



Api 



-A(s c )slu'(s c ), 



(28) 



e.g. A/iB ~ -1/Sc for a LJ potential. This asymp- 
totic behavior also holds for the more complicated cases 
discussed below. 

Mixtures. Many experimental relevant systems have 
mixed potentials such as the KA model for binary col- 
loidal mixtures sketched above. In general the interaction 
potential U a b(r) = u a b{s) between beads of two species 
a and b takes different energy parameters which causes 
different weights at the cutoff depending on which parti- 
cles interact. The impulsive corrections of such mixtures 
are readily obtained by linear superposition of Eq. (|25|) 
for different contributions (a,&). Let c a — p a /p denote 
the mole fraction of species a, cr a h the interaction range 
between a bead of type a and a bead of type b and g a b(s) 
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the respective radial pair distribution function. The im- 
pulsive correction to the Born-Lame coefficient thus be- 
comes 



f 



A ^B = -r 



f 



-p 2 A{s c )sl 



2d(d+l) 
^2^2c a c b ai b u' ab (s c ) g ab (s c 



(29) 



where we have used that for all types of interaction we 
have the same reduced cutoff s c . 

Let us now assume a mixture described by the gen- 
eralized LJ potential u a b(s) — e a b(l/s 12 — 1/s 6 ) with 
s = r i 'a a b. A reference energy e re f and a reference inter- 
action range cr re f may arbitrarily be defined using, say, 
the interaction of two beads of type a = b = f, i.e. 
e re f = en and a rc f = <j\\. Defining the dimensionless 
ratios w a b = £ab/^ef and v ab = {a ab /cr rc{ ) d we may thus 
rewrite the general Eq. ([2U1) as 

247r d / 2 

A W3 = - —j , T^FTTT^ P Xef e ref /lj(s c ) 



d(d + 2)r(d/2)' 
x c a CbVabW ab g a b(s c 



(30) 



Since g a b(s c ) 1 for large s c , the function /lj(s c ) de- 
termines the scaling as already stated above, Eq. (f28l) . 

Continuous polydispersity. We turn now to systems 
with a continuous polydispersity as in the second model 
investigated numerically below. Let us assume that each 
bead is characerized by a bead diameter D which is dis- 
tributed according to a well-defined normalized distri- 
bution d with t = D/o~ re f being a reduced bead di- 
ameter with respect to some reference length er re f. To 
be specific we shall assume a generalized LJ potential, 
Eq. (|21[) . where the interaction range aw and the energy 
scale ttt 1 between two beads are uniquely specified by the 
two reduced diameters t and t'. Defining wtt< = f-w Arcf, 
vtt' = (o~tt'/cret) an d using the radial pair distribution 
function gtt'(s) for two beads of reduced diameter t and 
t', the double-sum in Eq. (|30]l can be rewritten as the 
double-integral 

247^/2 

A ^B = - ,,, ; oxw./ox l erref e ref /lj(s c ) 



d(d + 2)r(d/2)' 
x dt dt' c t ct>v tt >wtt' gw{s c ). 



(31) 



In order to determine A/ib from Eq. (|31j) one needs to 
prescribe the laws for c t , a u > and e tt >. In the large-s c 
limit the double-integral becomes in any case constant, 
i.e. we have again Ap. B ~ -/lj(s c ) ~ -l/ s c~ d - 



III. COMPUTATIONAL ISSUES 

To illustrate the above predictions we present compu- 
tational data using two extremely well studied models of 
colloidal liquids at high temperatures which are described 
in detail elsewhere 123, \2M : 



Sc/so 


e/3 


Pp/p 


Kp/p 


Kp/p 


Gp/p 


Ap B p/p 


Kp/p 


0.9 


0.162 


4.61 


-19.2 


15.0 


-17.06 


17.08 


14.5 


1.0 


0.329 


5.39 


19.0 


19.0 


0.05 


0.03 


18.8 


1.1 


0.103 


4.90 


20.9 


17.2 


1.86 


-1.86 


17.2 


1.5 


-1.24 


3.13 


14.1 


13.2 


0.44 


-0.43 


13.3 


2.0 


-1.69 


2.71 


13.0 


12.4 


0.28 


-0.30 


12.3 


2.5 


-1.83 


2.55 


12.3 


12.1 


0.09 


-0.10 


12.0 


3.0 


-1.89 


2.50 


12.1 


12.0 


0.05 


-0.05 


11.9 


3.5 


-1.91 


2.47 


12.1 


12.0 


0.03 


-0.03 


12.1 


4.0 


-1.92 


2.46 


11.9 


11.8 


0.03 


-0.02 


11.9 



TABLE I: Various properties for polydisperse LJ beads at 
temperature T = 1 and density p ~ 0.72 vs. the reduced 
cutoff distance s c /so with sq = 2 1 / 6 being the minimum of 
the potential: energy per bead e, total pressure P — Pid + Pex, 
uncorrected compression modulus K, corrected compression 
modulus K — K+2ApB, bare shear modulus G and impulsive 
correction ApB obtained from the histogram fe(s) at s — s c . 
The corrected shear modulus G — G + AG vanishes as it 
should. The last column refers to the compression modulus 
K obtained using Eq. (|32l) for isobaric ensembles kept at the 
same pressure P (third column). 



• The already mentioned KA model [24[ for binary 
mixtures of L J beads in d = 3 has been investigated 
by means of Langevin MD simulation 0] imposing 
a temperature T = 0.8 for n = ?ia + «-b = 6912 
beads per simulation box, a total density p = 1.0 
and molar fractions c a = n^/n = 0.8 and q, = 
n-s/n = 0.2 for both types of beads A and B. As 
in Ref. [24[ we set ctaa = l.Oer, ctbb = 0.88cr and 
(tab = 0.8cr for the interaction range and eaa = 
l.Oe, ejjB = 0.5e and 6ab = 1.5e for the LJ energy 
scales. Only data for the usual cutoff s c — 2.5 is 
presented. 

• Using Monte Carlo (MC) simulations p], [|[ we have 
computed in d = 2 dimensions a specific case of 
the generalized LJ potential, Eq. (|2"Tj) . where all 
interaction energies are identical, e; = e, and the 
interaction range is set by the arithmetic mean 
ai = (Di + Dj)/2 of the diameters Di and Dj of 
the interacting beads. Following Ref. (2f| the bead 
diameters are uniformly distributed between 0.8a 
and 1.2cr. For the examples reported here we have 
used a temperature T = 1.0, n = 10000 beads per 
box and a density p w 0.72. 

We use LJ units throughout this work and Boltzmann's 
constant fee is set to unity. For the indicated parameter 
choices both systems correspond to isotropic liquids. The 
Table summarizes various properties for polydisperse LJ 
beads for different reduced cutoff distances s c /so. Con- 
sidering thermodynamic properties per particle (rather 
than per volume), we have made the data dimensionless 
by rescaling with the inverse temperature /3 and the den- 
sity p. 
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FIG. 1: Weighted radial pair distribution function fe(s)/3/p 
with s = r/ai being the reduced distance between two beads 
i and j. Main panel: KA mixtures in d — 3 (bold line) and 
polydisperse LJ beads in d = 2 (open symbols) for large re- 
duced distances s > so where the potential is attractive. The 
filled sphere corresponds to the shear modulus G computed 
using Eq. {18} for the KA system not taken into account the 
impulsive correction. Inset: Polydisperse LJ beads for s < so 
where fe(s) becomes strongly negative. 



IV. COMPUTATIONAL RESULTS 



A. Weighted pair distribution function /i2(s) 



The weighted pair distribution function fi2(s), 
Eq. (fT6|) . is presented in Fig. [1] Several cutoff distances 
s c are given for the polydisperse L J model, but for clarity 
only for distances s < s c . For the KA model only one 
cutoff is given, but this also for s > s c . Note that albeit 
different s c for each model correspond strictly speaking to 
different state points — as better seen from the energies 
per bead e or total pressures P indicated in the Table — 
the histograms vary only weakly with s c . Strong differ- 
ences become only apparent for very small s c as shown 
for s c = 0.9so in the inset. One can thus use the his- 
togram obtained for one s c to anticipate the impulsive 
correction for a different cutoff. Note that for large dis- 
tances corresponding to an attractive interaction we have 
h2(s) > (main panel). Obviously, h,2(s) vanishes at the 
minimum of the potential s = sq and for very large dis- 
tances s. Since g(s) ~ 1 in the latter limit, the histogram 
/i2(s) is given (up to a known prefactor) by s d+1 u' s (s). As 
one expects the decay is faster for the d — 2 data than 
for the KA mixtures in d — 3, since the phase volume 
at the cutoff is larger for the latter systems. Since all 
histograms are rather smooth, one may simply set s = s c 
for obtaining A/ib from /i2(s) instead of properly taking 
the limit s — >• s~Z . 



B. Compression modulus K 

The compression modulus K may be obtained from 
Eq. (|17[) or, equivalently, using the Rowlinson formula 
given elsewhere our systems are highly 

incompressible, i.e. the compression modulus K is large 
as usual in condensed matter systems, and it is thus dif- 
ficult to demonstrate the small correction predicted by 
Eq. (fT9|) . For the KA model we obtain, e.g. AK(3/p « 
(5/4) x 0.69 —1.2 which compared to the uncorrected 
estimate Kfi/ p s=s 21.9 is not very impressive. 

More importantly, it is not easy to obtain an inde- 
pendent and precise K-value for canonical ensembles of 
mixtures and polydisperse systems using, e.g., the total 
particle structure factor [13, lift- For polydisperse LJ 
beads we have thus computed K directly from the vol- 
ume fluctuations SV in the isobaric ensemble 111 



K = k B T 



(V) 
(6 2 V) 



(32) 



where we impose the same (mean) pressure P as for the 
corresponding canonical ensemble. As may be seen from 
the last column indicated in the Table, this yields similar 
values as the stress-fluctuation formula, Eq. (fTT|) . Un- 
fortunately, for larger cutoffs our error bars become too 
large to confirm the correction. The most striking exam- 
ple, where Eq. (TT9")) can be shown to work, is the case of 
the small cutoff s c = 0.9so : Using Eq. (TTT)) an impossible 
negative value K/3/p w —19.2 is obtained. As may be 
seen from the inset in Fig. [U one gets A^bP/p ~ 17.1 
from the weighted histogram /i2(s). Taking the correc- 
tion Eq. (TT9")) into account this yields Kf3/p w 15 which 
is similar to the value obtained using Eq. 



C. Shear modulus G 

Asymptotic limit for large sampling times. Since all 
our systems are liquids, the shear modulus G should of 
course vanish — at least in the thermodynamic limit for 
a sufficiently long sampling time. We have thus a clear 
reference and for this reason G is highly suitable to test 
our predictions. As can be seen from the solid sphere 
indicated in Fig. [T] for the KA mixtures with s c = 2.5 
we obtain Gj3/p ~ 0.65 if the impulsive correction for 
the Born term is not taken into account. As also shown 
by the figure, this deviation equals h2{s c )(3/ p w 0.69 as 
predicted. The same behavior is seen from Fig.[2]for poly- 
disperse LJ beads for a broad range of cutoff distances s c 
where the open squares refer to the uncorrected G and 
the filled squares to G obtained using Eq. ([2"(j]) . The solid 
lines indicated show Eq. (|3 1 1) . Focusing on the scaling for 
large s c we have set gtt>(s c ) = 1 in the double-integral 
which (under this assumption) is close to unity. 

Sampling time dependence. Figure [3] gives additional 
information for the shear modulus G(t) plotted as a func- 
tion of the number t of MC steps (MCS) for polydisperse 



7 




FIG. 2: Shear modulus G and impulsive correction — A/^b for 
polydisperse LJ beads vs. the reduced cutoff distance s c /so- 
The uncorrected shear modulus G (open symbols) has been 
obtained using the stress- fluctuation formula, Eq. (118[) . the 
correction term (stars) from the histogram /i2(s), Eq. (|16[1 . 
The solid lines indicate Eq. (|31|) where we have set <?tt'( s c) = 
1. Main panel: Linear representation showing that G = 
G+A(j,b (filled squares) vanishes as predicted, Eq. (|20[) . Inset: 
Double-logarithmic representation of the same data emphasiz- 
ing the asymptotic power-law decay for large s c as indicated 
by the bold dashed line. 
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FIG. 3: Shear modulus G for polydisperse LJ beads in d — 2 
for different s c as a function of the sampling time t given in 
units of MC Steps (MCS) of the local MC jumps used. The 
vertical axis is made dimensionless by means of a factor 13/ p. 
Filled symbols refer to the uncorrected G(t). The horizontal 
lines indicate — A/xb obtained from the histograms /i2(s) for 
three cutoffs as indicated in the Table. The dashed slope 
characterizes the decay of (the corrected) G(t) with time. 



LJ beads. Note that the sampling time is proportional 
to the number of configurations used for the averages. 
To smooth the data we have used gliding averages of 
length t over the total trajectories of length 10 7 MCS. 
For s c — l.Oso and s c = 4.0sq there was no need to add a 
correction while for s c = 0.9so where A/zb/3/p ~ 17.1 the 
uncorrected data is negative and cannot be represented. 



The filled symbols refer to the uncorrected shear modu- 
lus G(t) for s c = l.lso, s c = 1.5so and s c = 2.0so which 
are seen to approach for large times the predicted cor- 
rection — A/ib taken from the Table (horizontal lines). If 
corrected, all data sets vanish properly with time. 

Interestingly, neither /xb nor P ex do (essentially) de- 
pend on t while the fluctuation contribution — /xf(£) 
approaches (the corrected) hb — Pcx from below (not 
shown). The (corrected) shear modulus G(t) thus de- 
creases monotonously with time. As can be seen from 
Fig. [31 the (corrected) G(t) decays roughly as the power- 
law slope —1 indicated by the dashed line. (Note that 
the noise becomes too large for G(t)(3/p < 0.1.) Exactly 
the same behavior has been observed for the KA model 
in d = 3 (not shown). Apparently, G(t) decays quite 
generally inversely as the mean-square displacement h(t) 
of the beads in the free-diffusion limit, hit) ~ t. We re- 
mind that the same scaling G(t) ~ l/h(t) has also been 
reported for a bead-spring polymer model without im- 
pulsive corrections (s c — sq) 



V. CONCLUSION 

Summary. It has been emphasized in this study that 
an impulsive correction to the Born contributions Cg^ 7<5 
of the elastic moduli must arise if the interaction po- 
tential is truncated and shifted, Eq. ([3]), with a non- 
vanishing first derivative at the cutoff. To test our the- 
oretical predictions we have computed the elastic mod- 
uli of isotropic liquids in d = 3 and d — 2 dimensions. 
Since for these systems the shear modulus G must van- 
ish by construction, this allows a precise numerical ver- 
ification for different reduced cutoff distances s c . It has 
been shown how the impulsive correction for mixtures 
and polydisperse systems may be obtained from the read- 
ily computed weighted histogram /i2(s) which scales as 
~ s d+1 u'(s) for large s. As one expects, the cutoff 
effect vanishes if s c is large, Eq. ([28]) . or set to a mini- 
mum of the potential. It becomes more important with 
increasing dimension. 

Comment. Incidentally, it should be noted that the 
stress-fluctuation formula G = (1b + A*f ~ -fox and vari- 
ous other relations used in this work for liquid systems 
were originally derived for solids assuming well-defined 
reference positions and displacement fields @-H, [HI EH ■ 
It can be shown, however, that these assumptions can 
be relaxed and especially Eq. (fT5)) holds quite generally 
for isotropic systems. The aim of the present paper was 
to show numerically that the stress-fluctuation formal- 
ism yields the right value (G = 0) once the impulsive 
correction has been taken into account. 

Outlook. We are currently using the approach pre- 
sented here to characterize as a function of temperature, 
imposed pressure and sampling time the glass transition 
of the two models presented here. The generalization of 
our results to 

• other elastic moduli in anisotropic systems using 
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the more general impulsive correction Eq. (jlip . 
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